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We resolve discrepancies in previous analyses of the flow of collisionless dark matter particles in 
the Sun's gravitational field. We determine the phase-space distribution of the flow both numeri- 
cally, tracing particle trajectories back in time, and analytically, providing a simple correct relation 
between the velocity of particles at infinity and at the Earth. We use our results to produce sky 
maps of the distribution of arrival directions of dark matter particles on Earth at various times of the 
year. We assume various Maxwellian velocity distributions at infinity describing the standard dark 
halo and streams of dark matter. We illustrate the formation of a ring, analogous to the Einstein 
ring, when the Earth is directly downstream of the Sun. 



I. INTRODUCTION 

The motion of galaxies in the Universe cannot be explained by the gravitational pull of visible matter. Invisible 
dark matter (DM) is introduced to account for galactic rotation velocities and the gravitational interaction between 
galaxies. The spinning motion of the spiral disk of our Milky Way Galaxy is too fast to be explained merely by the 
gravity of its visible stars and gases, so it is assumed that the disk is surrounded by a much larger halo that contains 
a few stars but mostly unseen DM. A galaxy like our Milky Way is about 10% visible matter (stars and gas) and 90% 
unseen DM. 

The nature of DM is one of the unsolved problems of astrophysics and cosmology. Numerous theories have been 
proposed regarding the nature of the DM, from brown dwarf stars to MACHOs (Massive Compact Halo Objects) 
and theoretical subatomic particles such as WIMPs (Weakly Interacting Massive Particles) or axions. More than 20 
groups of physicists worldwide have been building devices to detect DM particles. 

Several of these devices aim to detect the elastic scattering of WIMPs off a laboratory target. Some of them also 
attempt to measure the direction from which the WIMPs reach the detector. These directional searches will most 
probably be the only observational way to explore the astrophysical properties of WIMPs near the Solar System. 
Current theoretical expectations for the local WIMP distribution range from a virialized Maxwellian velocity function 
(isothermal model) to a more or less complex network of streams of WIMPs. These streams arise either from the 
tidal disruption of satellite galaxies, like the Sagittarius dwarf, or from the process of gravitational collapse during 
the formation of the Milky Way Galaxy. 

For the planning and future interpretation of WIMP searches, it is important to identify and analyze the density and 
velocity distribution of DM particles in the Solar System. Damour and Krauss pj, Gould Q, Gould and Alam [J, and 
Lundberg and Edsjo 4] analyzed the transfer of DM particles from unbound to bound orbits caused by gravitational 
collisions with planets and by scattering with nuclei inside the planets or the Sun. Previously, Danby and Camm 
Danby and Bray @, Griest 0, and Sikivie and Wick @ studied the distribution of particles on unbound orbits, but 
their conclusions disagree with each other. All these papers assumed an isotropic Maxwellian distribution of velocities 
for the particles at infinity. 

To resolve the discrepancies in the latter studies just mentioned, and to prepare for an analysis of velocity distribu- 
tions more general than an isotropic Maxwellian, this paper is devoted to clarify the effect of the gravitational field 
of the Sun on the distribution of unbound orbits in our Solar System. We neglect the effects of the planets on the 
DM particle distribution. 

As an illustration of the problem we address, Fig. ^ depicts a flow of DM particles coming from infinity with a 
mean velocity V relative to the Sun and a velocity distribution function /(v) in the rest frame of the flow. The 
gravitational pull of the Sun deflects the particles, which speed up and get closer together as they approach the Sun. 
This so-called gravitational focusing effect changes the phase-space density of particles near the Sun. 

Our goal is to compute the phase-space density of dark matter particles /(r^.v^) as a function of the particle 
velocity v# at the position of Earth r^, starting from a given, but arbitrary, velocity distribution at infinity, i.e. far 
away from the Sun. 
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FIG. 1: Illustration of a flow of dark matter particles through the gravitational field of the Sun. 



The problem of the velocity distribution of a cloud of non-interacting particles in the gravitational field of a point 
mass was discussed by Danby and Camm [j| many years ago. They gave an expression for the velocity distribution 
function at any point in the Sun's rest frame, but confined themselves to points along the axis of symmetry. Few years 
later, the same problem was considered again by Danby and Bray 6j , where they studied the gravitational enhancement 
of the density of a non zero temperature interstellar matter near a star by applying the velocity distribution function 
of Danby and Camm 5]. In their paper, they discussed the problem away from the axis of symmetry. In 1988 
Griest studied the effect of the Sun's gravity on the distribution and detection of DM near the Earth in the case 
of the isothermal model. Griest |2| considered the annual modulation for the signal of DM detection and obtained 
his velocity distribution function from Danby and Camm |5j], correcting however some errors and transforming it to 
the Earth's rest frame. He found that the inclusion of bound particles in calculating the distribution function has 
negligible effect on direct detection of DM. In 2002 Sikivie and Wick Q analyzed the Sun's gravitational field for 
a flow with zero velocity dispersion through the Solar System. They gave expressions for the density and velocity 
distribution functions that they state were different from the analogous expressions in Danby and Camm [f| and 
Griest pj ■ In this paper we show that Danby and Camm's expression is correct provided their unspecified parameter 
/3 is set to -1, that Danby and Bray's expression need corrections, that Griest's expression is correct after fixing a typo 
evident by comparing other equations in his paper, and that Sikivie and Wick's expression is correct. All expressions 
but Danby and Bray's match our analytic formulas and numerical results. 

This paper is organized as follows. In Sec.[n] we give a general discussion of the distribution function for a flow of 
DM and its Boltzmann equation. In Sec. lIIII we describe a numerical backward-in-time method to determine the DM 
particle distribution function based on an adjustable time-step and a predictor-corrector iteration. Using this method, 
in Sec.|V]we show sky maps of the distribution of arrival directions of the particles on Earth. In Sec. II VI an analytical 
expression for the particle velocity at infinity Vqo is derived. When used to compute the particle distribution function, 
it gives the same results as in the numerical backward-in-time method. In Sec. lVIl we compare our results to the four 
analytical expressions for unbound orbits mentioned above. Finally, we summarize our results in Sec. IVIII 



In a flow of DM, particles move with a variety of velocities. Although it is hard to say what velocity an individual 
particle may have, it is possible to use a statistical approach to characterize the velocity attributes. The mathematical 
description of this statistical approach is called the phase-space distribution function f(r,v,t). 

For a system of N identical classical particles, the distribution function /(r, v, t) is given by 



where dN is the number of particles in the volume element d 3 x centered at r whose velocities fall in the velocity 
element d 3 v centered at v at time t. From the distribution function one can compute the macroscopic variables 
for an ensemble of particles. For example, the density p(r, t), the number of particles per unit volume, is given by 
p(r, t) = 6N/SV, where SN is the total number of particles in the volume element SV at r. In the limit SV — > d 3 x, 
SN is the integral of dN over all velocities. Hence 



II. THE DISTRIBUTION FUNCTION 



dN = f(r,v,t)d 3 xd 3 v, 



(1) 




(2) 



and 




(3) 
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An equation that governs the evolution of the distribution function / is the Boltzmann equation 

df df df fdf\ 



dt dxi dvi \dt 

Here Xi, Vi, and aj are the Cartesian components of the position, velocity, and acceleration vectors, respectively, and 
a summation is implied over repeated indices. The collision term on the right hand side receives contributions from 
particle collisions. However, cold dark matter is collisionless, thus 

For a stationary distribution function / we have df /dt = 0. In this case, the Boltzmann equation reduces to 

Vl d^ + ai d^ = °- (6) 

It follows that / is constant along trajectories, that is if r = r(t) and v = v(£) represent a particle trajectory, 
f(r(t), v(i)) = const, independent of t. This is Liouville's theorem. 

To find the value of f(jc E , v E ) at the position of the Earth, we determine the velocity Voo it had when it was far 
away from the Sun. Then we use Liouville's theorem to equate f(r E ,v E ) to the distribution function at infinity 

/oo(Voo)> 

f(r E ,v E ) = /oo(voo). (7) 

This is the basis of our analysis. 

It must be paid attention to the reference frame in which /oo(voo) is given. In our discussion, we assume it is in 
the frame of the Sun, namely that the velocity is measured relative to the Sun. For example, in the particle flow 
illustrated in Fig.^ which approaches the Sun with mean relative velocity V and with rest-frame distribution /(v), 
a Galilean transformation gives 

joo(Voo) =/(Voo- V), (8) 

In particular, for a Maxwellian distribution, the distribution function in the rest frame of the flow reads 

where is the uniform particle density at infinity and a is the velocity dispersion. In the Sun's frame, it becomes 

t i \ Poo f |voo-V| 2 \ 

- (^371 CX P ■ (10) 



III. NUMERICAL METHOD 



A backward-in-time method is used to compute the distribution function / for a flow of DM particles at the Earth, 
given a velocity distribution /oo(voo) far away from the Sun. 

To find the value of /(r^,VE) at the position of the Earth, we place a particle at the point (i\e, vg), and follow 
its trajectory numerically backward in time until it is far away from the Sun, at position with velocity (see 
Fig. 0). Since / is constant along every trajectory (Liouville's theorem), we have 

/(r B ,v £ ) = /(r 00 ,v 00 ). (11) 

We set the latter equal to the given velocity distribution /oo(Voo)- We choose jr^! =90 AU. 

The initial position ro of each simulated trajectory is taken to coincide with the position of the Earth t e , 

r Q =r E - (12) 
The initial velocity vq is set to the opposite of the arrival velocity of the particle on Earth, 



v = -v B , 



(13) 
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FIG. 2: The figure shows the idea of the backward-in-time method to compute the distribution function /. 



so that the particle follows its trajectory backward in time. 

For advancing particles along the trajectory, we use a modified version of the algorithm described in [9j. Newton's 
theory of gravity is used to produce a simulation of the motion of particles in the gravitational field of the Sun. In 
the simulation, the trajectories of DM particles in the Solar System are traced as they are deflected by the Sun's 
gravitational field. The accuracy of the numerical calculation is maintained by two techniques: an adjustable time-step 
for advancing the orbit, and a predictor-corrector iteration for improving the accuracy of each step. 

As is well-known, the acceleration of gravity on the particles due to the Sun (of mass M and located at the origin 
of the coordinate system) is given by 

a r = -r, 14 

where G is Newton's gravitational constant. 

To move the particle along the trajectory, the calculation uses a number of steps, up to 50,000, with a varying 
time-step At. If the acceleration were constant, equal to a, say, the velocity of the particle would change in a time At 
by Av = aAi, and its position would change by Ar = vAt with a velocity v equal to the average of the velocity 
at the beginning of the step and the velocity v.; + i at the end. Since the particle moves in a non-uniform gravitational 
field, there is an inevitable change in the acceleration. But if the change in position is sufficiently small, the error we 
make by assuming a constant acceleration can be made as small as we wish. 

Using the approximation that the position changes by the average velocity during the time-step and the velocity 
changes by the average acceleration, we have for the new velocity and position 

r l+ i = Yi + -fa + v i+ i)At, (15) 

v l+ i = Vj + i(a.; + a l+1 )At. (16) 

where a^ = a(r^), see Eq. (|14|) . This is a non-linear system of six equations in six unknowns, rj+i and Vi+j.. 

To compute v^+i from Eq. I|16|) we need to know the value of r^ + i that enters in a^+i, but to find r^+i we need to 
know v i+1 in Eq. i|15|) . To break the cycle, we use an iterative method, i.e. a series r^f-p v i+\ of successively more 
accurate approximations to the position and velocity at the end of the step (here k is the iteration number) . We start 
the iteration by setting 

r|+\ = r.+v.Ai+ia^Ai) 2 , (17) 

vg?! = v j: + a 2 Ai. (18) 

We refine the estimate using Eqs. I|15fl and (|16|) with and vji\ in the right hand side, and r^+i and vS.^ in 
the left hand side. Our convergence condition is 



(fe+i) _ ( fc ) 

V i+1 V i+1 



< e P \aiAt\ , (19) 



5 




where e p is a dimensionless number. After some trials, we found that e p — 10~ 3 was a good compromise between 
accuracy and efficiency. 

The time-step is adjusted as follows. At every step i, we compute the acceleration a i+1 at the new position using 
a trial time-step, and then compare the change in acceleration a,; + i — a; during the step with the acceleration aj at 
its beginning. If 

|a i+ i - aj| > e u |aj| , (20) 

where eu is a small positive number, we deem At too large. In this case, we divide At by two and restart the current 
step. Else, if 

|a s :+i - a*| < e ts W , (21) 

where e ts is another small positive number, we deem At too small. In this case, we double At and restart the current 
step. This is repeated until the time-step At is acceptable. After some experimentation, we found that e t ; = e ts = 1CP 3 
is a good choice. 



IV. ANALYTICAL METHOD 



As an alternative to the numerical method, and as a cross check of the calculation, the velocity at infinity 
appearing in Eq. (J7J can also be computed analytically. 

The expression for Vqo can be obtained using the conservation of the Laplace-Runge-Lenz vector following Sikivie 
and Wick We give here a shorter derivation that will lead to a simpler formula for v^. 

The Laplace-Runge-Lenz vector per unit mass is 

A = v x (r x v) - GMv = r±v 2 r ± - GMr. (22) 

Here rj_ = r±/r±, and rj_ is the projection of r perpendicular to v (see Fig.|3Jl- 
Conservation of the Laplace-Runge-Lenz vector A implies 

r±v 2 t±-GMr = bv 2 00 h + GMv 00 , (23) 

where b is the impact parameter and b is a unit vector perpendicular to Voo in the same plane as v and r (see Fig.|2J. 
Conservation of angular momentum in the form 

bvoo = r x v (24) 

allows us to eliminate b from Eq. I|23() . Dividing the result by r and using the identity r = +r • vv, the conservation 
of A equation becomes 

rx ( 2 GM\ _ GM. _ r± c GM . 

— [v r ± r-vv= — OToob-l v^. (25) 

r \ r J r r r 
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FIG. 4: Rotation of the basis vectors as a particle moves in the Sun's gravitational field. 

Eq. 1|25[) is a relation between the unit vector bases (fx, v) and (b, Vqo) (see Fig.0J). Introducing the rotation angle 
a between the two bases, we write 



b = cos ar j_ — sin av, 
Vqo = sinaf^ + cosav, 



Inserting Eqs. (|26[) and (|27|) into Eq. I|25|) . we have, in the (fj_, v) basis, 



r x GM . 

— Woo cos a H sin a 

r r 

GM rx 

cos a H vVac sin a 

r r 



GNI 



GM 



Solving this linear system for cos a and sin a, we find 



cos a 



sin a 



1 

D 



;2 _GMX _,GM ] r v 
r ) \ r 



1 r_L GM 
D r r 



GM 



ufooV ■ r 



where 



D = \ v — 



GM 



w M v • r v - 



GM 



- OTooV • r 



Energy conservation in the form 



2GM 



(26) 
(27) 



(28) 
(29) 

(30) 
(31) 

(32) 
(33) 



has been used in deriving the above expression for D. Notice in passing that D = \A\ 2 /r 2 . 

Substituting the above expressions for cos a and sin a into Eq. I|27[). multiplying by v^, and canceling common 
factors, we finally obtain the equation 



*4 V + Voo(GM/r)r - i; M v(v ■ f) 
^ + (GM/r)-v 00 (vf) 



(34) 



where Voo is given in Eq. I|33|) . As a simple check, when M = 0, Eq. (|34|l correctly gives — v. In our application, 
v and r are the particle velocity and position when it reaches the Earth, namely v = ve and r = ye- 
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FIG. 5: The ecliptic longitude of the Sun A Q changes from 0° to 360° as the Earth orbits the Sun. 

The expression for v M in Eq. (|34[) depends on the Keplerian character of the potential, but is independent of the 
velocity distribution function at infinity. This means that it can be used not only for a Maxwellian distribution but 
also for more general functions. 

When we used Eq. (f 34|) in the velocity distribution function, Eq. 1)10(1. we achieved an excellent agreement with 
our numerical backward-in-time method described in the previous section. Later, in Sec. IVII we will show that the 
expression for in Eq. (|34|l matches formulas used by Griest Q and Sikivie and Wick |J 

V. PHASE-SPACE DISTRIBUTION 

In this section, we present the velocity distribution of the DM particles as they arrive at the Earth. We place the 
Earth at four different locations during the year, and plot the flux of WIMPs as a function of their arrival direction 
n at fixed arrival speed ve- 

A. Positions of the Earth 

We specify the position of the Earth at different times of the year by means of the ecliptic longitude of the Sun A Q , 
which varies from 0° to 360° during the course of one year. Fig. |S] shows how A Q changes as the Earth orbits the Sun. 
At the Vernal (or Spring) Equinox (w March 21), A = 0°; at the Summer Solstice (« June 20), A Q = 90°; at the 
Autumnal Equinox (w September 21), A Q = 180°; and at the Winter Solstice (« December 20), A = 270°. Tabled 
summarizes the relation between the ecliptic coordinates of the Sun, the beginning of the four astronomical seasons, 
and their approximate calendar dates. 



TABLE I: The table relates the position of the Sun in the sky in ecliptic coordinates (Xq,/3q) to the beginning of the four 
astronomical seasons and their approximate dates. 



Sun's ecliptic coordinates (\q,[3q) 


Season 


Date 


(0°,0°) 


Vernal Equinox 


w March 21 


(90°, 0°) 


Summer Solstice 


w June 20 


(180°, 0°) 


Autumnal Equinox 


~ September 21 


(270°, 0°) 


Winter Solstice 


~ December 20 
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In the approximation of a circular orbit, one has the following relation between A© and the time of year t: 

3fif)° 

\ Q c^-(t~t VE ). (35) 
lyr 

Here t is the time during the year, which e.g. runs from at New Year's midnight to lyr ~ 365.2425 days at the end 
of the year. The time iyE is the time of the Vernal Equinox, i.e. iyE — 79.25 days from New Year's midnight to the 
midnight of March 21. 

We fix the origin of our coordinate system at the position of the Sun (see Fig.|3Jl. The x axis points in the direction 
of the Earth's position at the Autumnal Equinox, the y axis in the direction of the Earth's position at the Winter 
Solstice, and the z axis is perpendicular to the Earth's orbit. As seen from the Earth, the x axis points toward the 
position of the Sun at the Vernal Equinox, the y axis toward the position of the Sun at the Summer Solstice, and the 
z axis toward the North Pole of the ecliptic. 

The position vector of the Earth is given by (for a circular orbit) 

ye = — (cos A Q x + sin A Q y) AU. (36) 
Here 1 AU is the average distance between the Earth and the Sun. 



B. Velocities of the DM particles 



To ensure that only unbound orbits are considered, we constrain the particle speed ve to be equal to or greater 
than the escape speed v esc at the Earth's position, 



VE > Vcsc = \ • (37) 

V r E 

We specify the arrival direction of DM particles using a unit vector n pointing in the direction of arrival. We 
represent the sky by a sphere centered at the Earth (celestial sphere, see Fig. EJ. We use the ecliptic coordinate 
system in which the two coordinates of a point P on the celestial sphere are: (1) the ecliptic latitude 0, which is the 
angular distance from the ecliptic to P and varies from +90° at the Ecliptic North Pole to —90° at the Ecliptic South 
Pole, and (2) the ecliptic longitude A, which is the angular distance along the ecliptic from the Vernal equinox to the 
great circle through P and is measured eastwards in degrees from 0° to 360°. 

In ecliptic coordinates, 

n(A, (3) — cos P cos Ax + cos (3 sin Ay + sin /3z. (38) 

As can be seen from Fig. [5] since the unit vector n points in the direction from which the DM particle is coming, 
its velocity points in the opposite direction. Thus the velocity of the particle at the Earth is 

v E = —v E r\ = — ws(cos/3cos Ax + cos (3 sin Ay + sin/3z), (39) 

In the backward-in-time method, the initial velocity of the DM particle used to start the simulation is opposite to 
the actual velocity of the particle at the Earth, that is vq = — v_e = +«Eti. 



C. Sky maps of the DM particle flux 



In this Section, we present sky maps that show the arrival distribution of the DM particles for three cases: a stream 
moving in the +y direction, a stream moving in the —x direction, and the Maxwellian distribution of the standard 
halo. 

We start by deriving a formula for the flux of DM particles that reach the Earth. The number of particles at the 
Earth position is 

dN = f(r E ,VE)d 3 xd 3 v. (40) 
For a flux of particles, we write the volume element as 

d 3 x = dAdl = dA v E dt, (41) 
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FIG. 6: Specification of the arrival direction of a dark matter particle by means of a unit vector n and ecliptic coordinates A 
and P on the celestial sphere. 



where dA is an element of area orthogonal to the particle flux, and dl = VEdt is the length traveled by the particles 
of speed v E in a time dt. We use spherical coordinates in velocity space, and write the velocity element as 

d 3 v = v E dv B dQ., (42) 

where dfl is the solid angle in velocity space covered by the arrival directions of the particles. Inserting Eqs. I|41() and 
()42J) into Eq. H40|l , the number of particles becomes 

dN = f(r E ,v E )v%dAdtdv E dQ. (43) 

We define the specific flux of DM particles reaching the Earth from a direction n = — as the flux of particles per 
unit area per unit solid angle per unit speed 

F ^ h ^ dAdZ Ed n = f ^ v *- (44) 

We produce sky maps of F(v E , n) by fixing the value of v E and varying n on a regular grid of 360 points in A and 
100 points in sin/3. A gray scale is used to represent the values of F(v E , &)/poc, darker gray levels corresponding to 
lower values of the specific flux. In all the maps, the observer is on the Earth and looking toward the +x axis (see 
Fig. 0), and the arrival speed of the DM particles is fixed at v E = 30 AU/yr. 

Fig. shows the specific flux for a stream of DM particles moving toward the +y direction with velocity V y = 10 
AU/yr and V x — V z — 0. The stream velocity distribution at infinity is assumed to be gaussian, Eq. 1)11)0. with 
velocity dispersion a — 1 AU/yr. This value of a is chosen so as to reduce the dispersion of the flow and have a clear 
and sharp view of the distribution. The map in Fig. 0(a) is for A0 = 0°, which occurs at the time of the Vernal 
Equinox. In this map, the Sun is represented by a white dot at (A,/3) = (0°,0°). The specific flux is concentrated in 
the circular spot centered at (A, (3) — (270°, 0°), which is in the —y direction, as it should. The map in Fig. [If b) is for 
A© = 270°, which occurs at the time of the Winter Solstice. In this map, the specific flux has spread out into a ring 
around the Sun. The origin of this ring can be understood by referring to Fig. In the case we are discussing, the 
Earth is located at the point on the left of the Sun where the two trajectories drawn cross each other. An observer 
on the Earth sees particles coming from two directions in the plane of the figure. One is the flux of particles passing 
on one side of the Sun, the other is the flux of particles passing on the other side. When the figure is rotated around 
the Earth-Sun axis, the two directions trace the ring in Fig. 0b). This ring is analogous to the Einstein ring in the 
gravitational lensing of light. 

As a second example, we consider a stream approaching the Solar System from the +x direction, with V x = — 10 
AU/yr and V y = V z = 0. We again assume a Maxwellian velocity distribution at infinity with dispersion a = 1 AU/yr. 
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FIG. 7: Specific flux of dark matter particles F(ve, n) for a stream moving toward the +y direction at (a) Vernal Equinox, (b) 
Winter Solstice. The gray scale on the right shows the values of F(ve, n)/poo- 



The corresponding particle specific flux is shown in Fig. [5] at four successive times between the Vernal Equinox (« 
March 21) and the Summer Solstice (« June 20). At the Vernal Equinox, Fig.[|Ja), the Sun is directly in front of the 
stream, at A© = 0°. The Earth is on the leeward side of the stream, i.e. directly behind the Sun (stream, Sun, and 
observer are lined up in this order) . In this case, the observer sees the distribution of particles as a ring around the 
Sun, as explained before for Fig.[7{b). Ten days later, Fig. IHIb), the Sun has shifted eastward to X Q ~ 9?9. Now the 
observer sees an incomplete ring, because the Earth has moved to a point where fewer particle trajectories cross each 
other. As the observer's location changes again (see Fig. [SJc) for Aq ~ 19?7, 20 days after the Vernal Equinox), the 
ring starts to disappear because fewer and fewer trajectories cross. Eventually, the observer sees a circular spot, like 
in Fig. 03d) at A© = 90°, three months after the Vernal Equinox. 

Finally, we give an example of arrival distributions under the assumption of a standard dark halo. In this model, the 
DM particles are on average at rest relative to the Galaxy, and their velocity distribution is Maxwellian with velocity 
dispersion a — vlsr/v^- Here ulsr is the speed of the Local Standard of Rest (LSR), which moves at ulsr = 220 
km/s relative to the Galactic rest frame in the direction of the Galactic rotation. The latter is (l,b) = (90°, 0°) 
in Galactic coordinates, and (A,/3) = (347?340, 59?574) in ecliptic coordinates. In Astronomical Units, a — 32.816 
AU/yr. The other parameter we need is the mean velocity of the flow V with respect to the Sun. Alternatively, we 
can specify its opposite —V, i.e. the velocity of the Sun with respect to the flow. In the standard halo, the DM 
particles are on average at rest in the Galactic rest frame, thus —V is the velocity of the Sun with respect to the 
Galactic rest frame. We write it as the sum of the velocity vs un -LSR of the Sun with respect to the LSR (called 
the solar motion) and the velocity vt,sr of the LSR with respect to the Galactic rest frame. We take the galactic 
components of the solar motion to be |l(j U = 10.00 ± 0.36 km/s (radially inwards), V = 5.25 ± 0.62 km/s (in the 
direction of Galactic rotation) and W = 7.17+ 0.38 km/s (vertically upwards). The central values have a magnitude 
wSun-LSR = 13.378 km/s = 2.822 AU/yr in the direction (l,b) = (27?70, 32?409) or (A,/3) = (248?35, 32?189). Putting 
things together, we have 



V 



-VLSR — Vsun-LSR- 



(45) 



Converting from the Galactic to the ecliptic coordinate system and then our (x,y,z) coordinate system, we obtain 
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FIG. 8: Sky maps showing the arrival distribution of 30 AU/yr dark matter particles moving in the —x with velocity dispersion 
1 AU/yr. (a) A Q = 0° (Vernal Equinox), (b) A Q ~ 9?9, (c) A Q = 19?7, and (d) A Q = 90° (Summer Solstice). 
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V = (-22.049,7.372,-41.521) AU/yr in the direction (A,/?) = (-18?487, 60?755). 

The arrival distribution of DM particles for ve — 30 AU/yr at four times of the year is shown in Fig. [HJ As before, 
the intensity of the gray levels represents the values of the particle specific flux F(uE,n)/p 00 , with whiter regions 
corresponding to a higher flux. The effect of gravitational focusing near the Sun is evident. In order to help the 
reader place the DM particle distribution on the sky, the location of the 1986 brightest stars (brighter than ~ 5.17 
visual magnitude) is indicated by black dots. The size of each dot is proportional to the star's magnitude. Notice the 
constellation of Orion at A w 85°, (3 w —30°, and the constellation of Ursa Major at the top left of the figure. An 
almost sinusoidal "band" of stars can be followed from the lower left corner through Orion up to the top center and 
symmetrically down: it is the Milky Way. 



VI. COMPARISON WITH PREVIOUS ANALYSES 



In this section, we compare our results with analytical expressions of Danby and Camm [5j, Danby and Bray |6J, 
Griest Q, and Sikivie and Wick 8]. We find Danby and Camm's expression correct after fixing their unspecified 
parameter /?, Danby and Bray's expression incorrect unless few signs are changed, Griest 's expression correct after 
fixing a typo, and Sikivie and Wick's expression correct. We checked that all distributions, with the exception of 
Danby and Bray's, match the results of our numerical and analytical methods shown in Figs. HOD and EH 



A. Danby and Camm 

In the work of Danby and Camm [j| , the velocity distribution function is determined at all points around the Sun. 
Cylindrical polar coordinates [w, ip, z) are used with the origin at the center of attraction, and the z axis directed 
towards the streaming cloud of particles (see Fig. llOfl . The coordinate w is the distance of the Earth from the flow axis, 
and if is the azimuthal angle. The corresponding velocity components are II, Z. Danby and Camm's distribution 
function at any point reads, in their notation, 



where 




c 2 + h + 2c H 1 p jyf^ — !—± , 46 

h + pr- 1 + pil'Z cos 9 + pil /2 Usin9 



h = U 2 + <P 2 + Z 2 - — . (47) 



Since II, and Z are the velocity components and p, = GM, then Ii is identical to our v 2 ^ given by Eq. (|33|) . 

For the parameters h and po, we have h = l/a and pa = p^. The parameter c in Eq. (|46|l can be obtained from 
the identity 

n 2 + $ 2 + (Z + c) 2 = |v- V| 2 , (48) 

which follows by comparing the exponents of the Maxwellian velocity distributions in Eqs. I|46|l and (JTUJ. We find 
c = V. 

For the other parameters in Eq. Q46f) . we use Fig. 1101 to define them in terms of our parameters and according to 
our notation. In Fig. 1101 Z is the component of the particle's velocity v along the z axis, II is the radial component 
of v (orthogonal to the z axis), and 9 is the angle between the +z axis and the position vector r. (The angle A is used 
by Griest and will be discussed later when comparing our calculation to his results.) Since the velocity of the Sun 
v s points in the +z direction, the angle 9 is also the angle between r and v s . On the windward side of the axis of 
symmetry (their z > 0), sin# = and cosd = +1. On the leeward side (their z < 0), sin 9 = and cos9 = —1. From 
the rotational symmetry of the problem, the velocity distribution is independent of the azimuthal angle <p. 

Using V = — v s = — z, we have 

cosfl = f • z = -f • V (49) 

and 



Z = v ■ z = —v • V. 



(50) 
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FIG. 9: Sky maps showing the arrival distribution of 30 AU/yr DM particles using standard halo parameters for a and V. (a) 
A = 0° (Vernal Equinox), (b) A Q = 90° (Summer Solstice), (c) A Q = 180° (Autumnal Equinox), and (d) A = 270° (Winter 
Solstice) . 
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Z<0 

leeward side 



z>0 

windward side 



FIG. 10: Relating Danby and Camm's and Griest's parameters to our parameters r, v, and V. 



Since II is the radial component of v, we can write it as 

n = ve CT , (51) 
where e ro is a radial unit vector orthogonal to the z axis. An expression for e ro valid off the flow axis is 



f VV 
f • VVI 



(52) 



The second equality follows from Eq. I|49|) and the fact that sinf? > for < 8 < tt. On the axis, sin 9 = and 
Eq. H52fl contains a division by zero. However, Danby and Camm's expression contains only the product e ro sin#, 
which is well-defined. Indeed, their expression contains the combination Zcos9 + II sin 9, which with the help of 
Eqs. (|5T|) and 1(521) becomes 



Z cos 9 + II sin 9 = — v ■ V cos ( 



v ■ f + cos 9 v • V 



sin ( 



According to our analysis and notation, Eq. 14rj|) therefore takes the form 

f _ ( F DC \ 

1 - (2^)3/2 exi \ 2a V' 



where 



i^c - ^ 2 + t& - 2V • 



^4v - 0v^{GM/r)r + Awv(y ■ f) 
^ + (GM/r) + fc(v.f) 



We used 71=^, and cos6>, Z, and Zcos6» + IIsm6> given in Eqs. (HHJ), JSOI), and (jSHJ - 

The quantity -Foe? can be put into the form 



F DC = V 2 + v 2 x - 2V • v 



oo — v oo 



(53) 



(54) 



(55) 



(56) 



with Vqo given by Eq. I|34(l provided the parameter (3 = — 1. 

Danby and Camm introduced [3 in Eq. (|46|l to deal with an ambiguity in the derivation of their formula. This 
ambiguity lead Danby and Camm [f| to confine their attention to the axis of symmetry instead of making full use of 
their expression. They solved their ambiguity by putting cos# = 1 and choosing (3 — +1 on the leeward side of the 
z axis and (3 = — 1 on the windward side. This choice gives the correct solution on the z axis, because on this axis 
sin# = and only the product (3 cos 9 appears in Eq. (|46|l . With (3 = — 1, the windward side has (3 cos 9 — — 1 and 
the leeward side (3cos9 = +1. These are the same values that Danby and Camm obtain on the axis. 

We conclude that Danby and Camm's expression is correct provided j3 = — 1 and care is taken in using Eq. (|53|l for 
Z cos9 + IIsin#. Once this is done, their expression gives the same results as our numerical and analytical methods 
shown in Figs. □ El and El 
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B. Danby and Bray 

Danby and Bray's distribution function, Eqs. (2) and (3) of their paper 0, reads, in their notation, 

/ = (27r)- 3 / 2 Pa h 3 expf^FuB^j , (57) 

where 

2 , t2 , - J 2 Z- J(p/r) cos 9 + JZ(Zcos9 + UsinO) 
F DB -c +J + 2c j2 + (M/r) + J(Zcos0 + nsin0) • (58) 

Danby and Bray's notation is the same as Danby and Camm's, except that J — I\ 1 and that the direction of the z 
axis in Fig. has been reversed. 

Danby and Bray's formula is incorrect, as we now show. In our notation we have [cfr. Eqs. 149l) . 15Ufl . and (|53|) ] 

cos0 = ?-V, (59) 

Z = v ■ V, (60) 

Zcos0 + nsin6» = v • f . (61) 

Hence, 

F0B - V +V -- 2V - v^ + iGMM+v^-r) ■ (62) 

This formula is incorrect because the v^v term in the numerator and the Uoo(v • f) term in the denominator have the 
wrong sign. To fix Danby and Bray's formula, one should change the sign of the J 2 Z term in the numerator and of 
the J (Z cos 9 + ITsin^) term in the denominator in Eq. iJSSJ. 



C. Griest 

Griest used spherical coordinates with the Sun at the origin and the positive z axis on the line from the Sun to 
the Earth. He considered a cloud of DM particles that, far away from the Sun, had a uniform density and an isotropic 
Maxwellian distribution of velocities. Griest's distribution function reads, in his notation, 

/ = ( 2 7r)- 3 / 2 ^ 3 exp(^F G ) 9 (J 2 ), (63) 

where 

2 2 J 2 Z + J(GM G / a(B ) cos 6 - JvZ cos A 

F G =v s +J +2v s — r . (64) 

J 2 - + (GM /a ffi ) — Jv cos A 

A typing mistake in Griest's paper was corrected in writing the above expression for Fq, namely cos A in the last 
term in the denominator was mistakenly written as cosO. This typing mistake can be seen by comparing Eqs. (2) 
and (5) of Griest's paper 0. 

The step function 9(J 2 ) in Eq. (gSJl incorporates the assumption that there are no particles in bound orbits. In our 
calculation, we always have 9 (J 2 ) = 1. 

In the expression of Fg, v s is the Sun's velocity with respect to the cloud of DM, and according to our notation 
v s = V. The parameter J 2 is identical to v 2 ^, which is given by Eq. Ij33(l . The parameters cos 9 and Z are the same 
as in Danby and Camm's paper [5J and related to our parameters by Eqs. (|49J) and (|50|) . Furthermore, the parameter 
v is the arrival speed of the DM particles on Earth with respect to the Sun, a§ is the distance between the Earth and 
the Sun, h is the inverse of the velocity dispersion, and po is the density of the particles at infinity. Converting to our 
notation, we have Mq = M, a§ = r, h = 1/ct, and po = poo. Finally, from Fig. HUI we can write 

cos A = v • f . (65) 

According to our analysis and notation, Griest's formula Eq. (|63|l takes the form 



16 



where 

Fg = V 2 + v 2 — 2V • ^ V + v ^ GM l r ^ ~ "°° v ( v • f ) (67) 

In writing the last expression, we used J = t>oo, and cos#, and cos A given in Eqs. I|49|). I|5U|I . and l|65|) . 
Eq. (|BT|) is equivalent to 

F G = V 2 +i4-2V- Voo = | Voo - V| 2 (68) 

provided Vqo is given by Eq. i|34|) . Thus Griest's distribution function is identical to our expression for the Maxwellian 
distribution function, Eq. (|10fl . and it gives the correct distribution as in Figs.[7||Hl and|5J 



D. Sikivie and Wick 



We also compared our results to the work of Sikivie and Wick Sikivie and Wick assumed that the velocity 
distribution in the galactic halo is isothermal. They gave a position-dependent phase-space distribution in the presence 
of the Sun which in their notation reads 



/(r,v) = 



■ exp 



'sw 



-v 2 G (r,v) 



©(*4(r,v))- 



(69) 



As in Griest's paper Q, the step function 9(w 2 c (r,v)) incorporates the assumption that there are no particles in 
bound orbits. In our calculation, we always have 0(w 2 c (r, v)) = 1. = is the density of the particles at infinity, 
and asw = y/2/3a. The factor v G (r, v) in the exponent of Eq. ljf)9"|) is given by || 



y c( r , v ) = ( v © + Voo(r,v)) 2 , 



(70) 



where v Q = —V is the Sun's velocity with respect to the mean velocity of the DM particles, and v 00 (r, v) is given by 
the following formula in Sikivie and Wick's paper: 



Voo(r,v) 



1 



a 2 ?; 2 + 1 2 



v[l 



av^r — av Q 



(r • v)] + raw 2 



1. . v 2 
-(r • v H 



(71) 



Here a = GM/v 2 ^, I 2 = r 2 v 2 — (r ■ v) 2 , and Voo is given in Eq. J22J). 

The distribution function, Eq. I|69|l . matches our DM particle distribution. Indeed, Sikivie and Wick's expression 
for Vqo, Eq. iJHJ, is analytically identical to ours, Eq. I|34|l . due to the cancellation of the common factor 



GM 



OT„oV • r, 



(72) 



between the numerator and the denominator of their expression. 



VII. CONCLUSIONS 



In this paper, we analyzed the distribution of a flow of unbound collisionless DM particles in the Solar System, 
both numerically and analytically. In particular, we focused on the velocity distribution at the location of the Earth. 
We used Liouville's theorem to relate the phase-space distribution at the Earth to the velocity distribution at infinity 
(SecHIJ. 

In the numerical method fSec. IHIfl . we traced the trajectories of the DM particles backward in time after they are 
deflected by the Sun's gravitational field. This numerical method is independent of the special form of the gravitational 
field and of the velocity distribution at infinity. The calculation used Newtonian gravity for advancing the trajectories 
until the particles return to infinity. The accuracy of the calculation was maintained by an adjustable time-step and 
a predictor-corrector iteration for improving the accuracy of each step along the trajectory. 

In the analytical method fSec. UVfl we obtained a formula, Eq. I|34|) . for the velocity of the DM particles at infinity. 
That formula is valid for motion in a Keplerian field, and is independent of the choice of velocity distribution at 
infinity. 
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We applied both the numerical and the analytic method to a selection of Maxwellian velocity functions at infinity 
representing streams of DM and the standard halo (Sec.0. Comparison of the numerical and analytical calculations 
gave the same results. 

For these velocity functions, we produced a number of sky maps showing the arrival distribution of DM particles 
at the Earth - actually the specific flux defined in Eq. (|44l) - at different times of the year (see Figs. and [8J. The 
maps also show the location of the Sun in the sky as the Earth moves around the Sun. The arrival distribution of DM 
particles under standard dark halo assumptions is displayed in Fig.0 which also shows the positions and magnitudes 
of the 1986 brightest stars in the sky. 

The effect of the gravitational field of the Sun on the distribution of DM particles is evident on these maps. For 
example, the gravitational deflection produces a ring in the arrival distribution when the Earth is directly on the other 
side of the Sun as seen from the stream (see Fig.[7Jb)). This ring is analogous to the Einstein ring in the gravitational 
lensing of light. 

Finally, comparing our results with previous analyses, we were able to resolve the issue of discrepant results in 
Danby and Camm Danby and Bray 0], Griest 0, and Sikivie and Wick 0. Danby and Camm's distribution 
function is correct on the axis of symmetry of the flow, and their ambiguity in the choice of the parameter (3 can be 
fixed by choosing (3 — —1. Danby and Bray's have an incorrect distribution with a couple of wrong signs. Griest's 
distribution is correct after fixing a typo evident from comparing various equations in his paper. Sikivie and Wick's 
distribution is also correct. We checked that Danby and Camm's, Griest's, and Sikivie and Wick's distributions match 
the results of our numerical and analytical expressions shown in Figs.[7||Hl and [5] 
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